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An approach is presented for theoretical calculations of the Debye- Waller factors in x-ray ab- 
sorption spectra. These factors are represented in terms of the cumulant expansion up to third 
order. They account respectively for the net thermal expansion cr'^-'(r), the mean-square relative 
displacements a^{T), and the asymmetry of the pair distribution function cr''^'(r). Similarly, we ob- 
tain Debye- Waller factors for x-ray and neutron scattering in terms of the mean-square vibrational 
amplitudes u^{T). Our method is based on density functional theory calculations of the dynamical 
matrix, together with an efficient Lanczos algorithm for projected phonon spectra within the quasi- 
harmonic approximation. Due to anharmonicity in the interatomic forces, the results are highly 
sensitive to variations in the equilibrium lattice constants, and hence to the choice of exchange- 
correlation potential. In order to treat this sensitivity, we introduce two prescriptions: one based on 
the local density approximation, and a second based on a modified generalized gradient approxima- 
tion. Illustrative results for the leading cumulants are presented for several materials and compared 
with experiment and with correlated Einstein and Debye models. We also obtain Born-von Karman 
parameters and corrections due to perpendicular vibrations. 
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I. INTRODUCTION 

Thermal vibrations and disorder in x-ray absorption 
spectra (XAS) give rise to Debye- Waller (DW) factors 
varying as exp[-VF(r)], where W{T) « 2Pa^{T) and 
cr'^{T) is the mean square relative displacement (MSRD) 
of a given multiple-scattering (MS) pathii These Debye- 
Waller factors damp the spectra with respect to increas- 
ing temperature T and wave number k (or energy), and 
account for the observation that the x-ray absorption fine 
structure (XAFS), "melts" with increasing temperature.— 
The XAFS DW factor is analogous to that for x-ray 
and neutron diffraction or the Mofibauer effect, where 
W{T) = {l/2)k^u^{T). The difference is that the XAFS 
DW factor refers to correlated averages over relative dis- 
placements, e.g., CT^ — {[{ur — uo) • R]^) for the MSRD, 
while that for x-ray and neutron diffraction refers to the 
mean-square displacements v?{T) = ((u • R.)^) of a given 
atom. Due to their exponential damping, accurate DW 
factors are crucial to a quantitative treatment of x-ray 
absorption spectra. Consequently, the lack of precise 
Debye- Waller factors has been one of the biggest limi- 
tations to accurate structure determinations (e.g., coor- 
dination number and interatomic distances) from XAFS 
experiment. 

Due to the difficulty of calculating the vibrational dis- 
tribution function from first principles, XAFS Debye- 
Waller factors have, heretofore, been fitted to experimen- 
tal data or estimated semi-empirically, e.g., from cor- 
related Einstein and Debye models.^''* However, these 
ad hoc approaches are unsatisfactory for several reasons. 
First, there are often many more DW factors in the MS 
path expansion than can be fit reliably. Second, semi- 
empirical models typically ignore anisotropic contribu- 
tions and hence do not capture the detailed structure of 



the phonon spectra. 

To address these problems, we introduce first principles 
procedures for calculations of the Debye- Waller factors in 
XAS and related spectra. Our approach is based primar- 
ily on density functional theory (DFT) calculations of 
the dynamical matrix, together with an efficient Lanc- 
zos algorithm for the projected phonon spectra.— »^ DFT 
calculations of crystallographic Debye- Waller factors and 
other thermodynamic quantities have been carried out 
previously using modern electronic structure codes;^!^ 
and our work here builds on these developments, with 
particular emphasis on applications to XAS. 

Due to intrinsic anharmonicity in the interatomic 
forces, the behavior of the DW factors is extremely sen- 
sitive to the equilibrium lattice constant a. For exam- 
ple, we find that varies approximately as a^^, where 
7 = —dlnuj/dhiV is the mean Griineisen parameter 
which is typically about 2 for fee metals, and uj refers 
to the mean phonon frequency. Consequently cr^ is also 
very sensitive to the choice of the exchange-correlation 
potential in the DFT, since a 1% error in lattice constant 
yields an error of 67 « 10% in . As a result, relatively 
small errors in the lattice constant predicted by the local 
density approximation (LDA) which tends to overbind, 
or the generalized gradient approximation (GGA) which 
tends to underbind, become greatly magnified^'^ in DW 
calculations. 

In order to treat this sensitivity we have developed two 
ad hoc prescriptions for ah initio calculations of Debye- 
Waller factors based on DFT calculations with I) the con- 
ventional LDA and II) a modified-GGA (termed hCGA) 
described below. For comparison we also present se- 
lected results with a conventional GGA, with the corre- 
lated Einstein and Debye models, and with an empirical 
model based on the Born-von Karman parameters ob- 
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tained from fits to phonon spectra. Detailed results are 
presented for a number of fee and diamond structures. 



II. FORMALISM 
A. Cumulants 

In this section we outline the formalism used in our 
approach. Physically, the DW factors in XAS arise from 
a thermal and configurational average of the XAS spec- 
tra over the pair (or MS path length) distribution 
bmction, where is the x-ray absorption coefficient 
in the absence of disorder. The effects of disorder and vi- 
brations are additive, but since the factors due to config- 
urational disorder are dependent on sample history and 
preparation, in this paper we focus only on the thermal 
contribution. The effect of the DW factors on the XAFS 
x(fc) is dominated by the average over the oscillatory be- 
havior of each path in the multiple-scattering (MS) path 
expansion xr{^) oc sin(2fci? -I- <&). If the disorder is not 
too large the average is conveniently expressed in terms 
of the cumulant expansion, l^ii 
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where r is the instantaneous bond length, i?o the equilib- 
rium length in the absence of vibrations, and ct*^"^(T) the 
n-th cumulant average. For multiple-scattering paths, 
this length refers to half the total MS path length. The 
dominant effect on XAFS amplitudes comes from the 
leading exponential decay factor W{T) « 2k^a^(T), 
while the imaginary terms in W{T) contributes to the 
XAFS phase. The leading such contribution is the 
thermal expansion which comes from the first cumulant 
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Thus, the mean bond length is (r) = f = Rq + a^^\T). 
The skew of the distribution, which is given by the third 
cumulant a^^^{T) contributes a negative phase shift, and 
hence the mean distance obtained in fits to XAFS exper- 
iment is typically shorter than that obtained from the 
first cumulant alone.-i As emphasized above, an accu- 
rate account of the effects of anharmonicity is key to a 
quantitative treatment of these DW factors over a broad 
range of temperatures. This is illustrated in Fig.[Tl which 
shows the strong variation in the mean phonon frequency 
D = u;/2iT vs small variations in lattice constant as calcu- 
lated using various models described below. The expres- 
sions for the higher cumulants in Eq. ([T]) simplify when 
expressed with respect to the mean, and are^ 
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FIG. 1: Mean frequency z? of the VDOS projected along the 
nearest neighbor single scattering path of Cu, obtained from 
the first Lanczos iteration. The vertical line indicates the ex- 
perimental lattice constant at 298 K while the horizontal line 
shows the Einstein frequency obtained from the experimen- 
tally determined DW factor. The Born-von Karman parame- 
ters for Cu at 298 K were taken from Ref. \li 



The thermal averages involved in the calculation of the 
cumulants can be expressed in terms of the projected vi- 
brational density of states (VDOS) pij(w).— i^^^^ For ex- 
ample, the MSRD for a given path R is given by the 
Debye integral 



liT) 



h 



, ff3huj\ , . , 
■ coth — — \pR (w) dw, 
ui V 2 / 



(5) 



where [ir is the reduced mass associated with the path, 
/3 = l/ksT, and pR [lj) is the vibrational density of states 
projected on R. In the following, the path index subscript 
R is suppressed unless needed for clarity. 

The first cumulant a^^^ is generally path-dependent 
and reflects the anharmonic behavior of a system. For 
monoatomic systems, this quantity is directly propor- 
tional to the net thermal expansion Aa = a{T) — ao, 
which can be obtained by minimizing the vibrational free 
energy F{a,T). Within the quasi-harmonic approxima- 
tion, F{a,T) is given by a sum over the internal energy 
E{a) and the vibrational free energy per atom 

F{a,T) =E{a) + 
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where T is the temperature, Pa{oj) is the total VDOS, 
and we have assumed cubic symmetry for simplicity. 

Furthermore, as pointed out by Fornasini et al.,^^ the 
values of the cumulants measured in XAFS experiments 
include two further corrections. First, perpendicular vi- 
brations lead to a small increase in the mean expansion 
observed in XAFS compared to that in x-ray crystallog- 
raphy: 
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FIG. 2: Total vibrational density of states of Cu at 49 K 
from DFT calculations using the LDA with our prescription 
I; the hGGA (see text) with our prescription II; and from 
experimen1ii^. 



We have shown (see Appendix) that cr^ and are closely 
related, and hence that cr^ can be estimated in terms of 
tj^ . Second, the position dependent XAFS amplitude fac- 
tors exp(— 2i?/A)/i?^ give rise to an effective radial dis- 
tribution function g{R) g{R) exp{—2R/ X)/ which 
shifts thermal expansion observed in XAFS by an addi- 
tional correction 
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This second correction is often included in XAFS analysis 
routines and has been taken into account in the experi- 
mental results presented here Note that the corrections 
in Eq. ([7]) and ([8]) are both of the same order of magni- 
tude and partially cancel. 



B. Lanczos Algorithm 

The VDOS pn (w) has often been approximated by 
means of Einstein and Debye models based on empiri- 
cal data. Although these models are quite useful, espe- 
cially for isotropic systems such as metals without highly 
directional bonds, their limitations are well knownJ^ii^ 
To overcome some of these limitations Poiarkova and 
proposed a method in which the VDOS is cal- 
culated from the imaginary part of the lattice dynamical 
Green's function 
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Here |0) is a Lanczos seed vector representing a nor- 
malized, mass- weighted initial displacement of the atoms 
along the multiple-scattering path R, and D is the dy- 
namical matrix of force constants 
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where Ujia is the a = {x, y, z} Cartesian displacement 
of atom j in unit cell I and Mj is the mass of atom j, 
and where E{a) is the internal energy of the system eval- 
uated at the lattice constant a{T). Thus, our approach 
takes into account the main effects of anharmonicity in 
terms of force constants that depend parametrically on 
the temperature. 

Efficient calculations of the lattice dynamical Green's 
function can be accomplished using a continued fraction 
representation, with parameters obtained with the itera- 
tive Lanczos algorithm. This yields a many-pole repre- 
sentation for the VDOS which is well suited for accurate 
spectral integrations. The first step in the Lanczos algo- 
rithm corresponds to the correlated Einstein model. 
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with an Einstein frequency uj given by 

L^2 = (0|D(i?)|0). (12) 

The frequency u corresponds to the rms average over the 
projected phonon spectra pito). However, the choice of 
the Einstein frequency is not unique, and the appropri- 
ate choice depends on the physical quantity being calcu- 
lated, as discussed in more detail below. Poiarkova et al. 
truncated the continued fraction at the second tier (i.e. 
second Lanczos iteration), which is usually adequate to 
converge the results to about 10%. Subsequently Krappe 
and Rossner— showed that at least six Lanczos iterations 
are required to achieve convergence to within 1%. Thus 
the Lanczos algorithm provides an efhcient and accurate 
procedure for calculating MS path-dependent DW factors 
from Eq. 

The main difficulty in implementing the Lanczos algo- 
rithm lies in obtaining an accurate model for the dynam- 
ical matrix (or force constants) D for a given system. Al- 
though semi-empirical estimates of interatomic force con- 
stants or Born-von Karman parameters are sometimes 
available, their temperature dependence limits their ac- 
curacy and generality. Similarly, simple models for the vi- 
brational distribution function (e.g., Einstein and Debye) 
generally ignore anisotropic behavior. One of our main 
aims in this paper is to develop a first principles approach 
that allows us to calculate the force constants for vari- 
ous systems using DFT. In addition, we have extended 
the Lanczos algorithm described above to several other 
cases by generalizing the Lanczos seed-state |0). This 
allows us to calculate several other quantities including 
the total vibrational density of states (VDOS), the vibra- 
tional free energy, thermal expansion, the mean square 
atomic displacements u^iT) in crystallographic Debye- 
Waller factorsj^i In addition, we calculate cr^, which 
yields the perpendicular motion contribution to the DW 
factor of Eq. and estimates of the third cumulant. 
Representative results for the VDOS calculated by this 
method are illustrated in Fig. [51 
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FIG. 3: Temperature dependence of the Debye- Waller factor 
for the nearest neighbor single scattering path in Cu. The 
experimental difference values^^ were shifted to match the 
LDA (I) results at K. 



FIG. 4: Temperature dependence of the Debye- Waller factor 
for the nearest neighbor single scattering path in Ge. The 
experimental difference values^S were shifted to match the 
LDA (I) results at K. 



Correlated Einstein Model 



III. DFT CALCULATIONS 



Although the cumulants other than the second are of- 
ten negligible for small anharmonicity, their calculation 
using the apparatus of anharmonic lattice dynamics is 
computationally demanding. On the other hand, it has 
been shown that these cumulants can be approximated to 
reasonable accuracy using a correlated anharmonic Ein- 
stein model for each MS p&th.^^ and this is the method 
adopted here. In this approach an Einstein model is 
constructed for each MS path keeping only cubic anhar- 
monicity, yielding the effective one-dimensional potential 



1 



(13) 



where x is the net stretch in a given bond. The Ein- 
stein frequency uje within the quasi-harmonic approxi- 
mation is then obtained from the relation Eq. (|A.4p . i.e., 
k — ko + Gk^x ^ Huj'^. This choice of Einstein frequency 
ensures that the high temperature behavior of from 
the Einstein model agrees with Eq. The construc- 
tion of this Einstein model from the dynamical matrix D 
along with explicit examples is given in the Appendix. 
The relations between the cumulants for the Einstein 
model can be used to obtain estimates for cr'-^-' and a'^^K 
For example, for the first cumulant 
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Note that this relation differs from that in Refs. 

in that it contains an extra multiplicative factor 77 = 

l/(aj~^)(D^, as discussed in the Appendix. 



A. Computational Strategy 

As noted above, one of the main aims of this paper is to 
calculate the force constants within the quasi-harmonic 
approximation using DFT and an appropriate choice of 
exchange-correlation functional. Due to the extreme sen- 
sitivity of the phonon spectra to the interatomic dis- 
tances, as discussed above, the most important parame- 
ters entering the calculation of the dynamical matrix are 
the lattice constant and the geometry of the system. A 
typical example of the effect of expansion is illustrated in 
Fig. [Jl which shows the variation of the first moment of 
the VDOS (i.e. the average frequency v) projected along 
the nearest-neighbor single-scattering path of Cu. For 
comparison Fig. [1] also shows v obtained with a model 
based on the Born- von Karman parameters at 298 K. As 
expected, when the system expands the vibrational fre- 
quencies are red-shifted due to the weakening of the inter- 
atomic interactions. From the common slopes in Fig. 1., 
we see that all of the functionals have similar Griineisen 
parameters 7 ~ 2.2 at the experimental lattice constant 
3.61 A, in accord with the experimental valued 2.0 ±0.2. 
Note that although at a given lattice constant the GGA 
functional always produces a stiffer model than LDA, i.e., 
with higher mean frequencies, the results at the equilib- 
rium GGA lattice constant tend to be softer than at the 
equilibrium LDA lattice constant. Moreover, when com- 
pared with the experimental value, the LDA and GGA 
functionals respectively underestimate and overestimate 
the mean frequency by about 5%. This translates into 
a 20-25% error in the DW factors calculated with these 
methods. This margin of error is too large to make the 
DW factors of significant value in quantitative EXAFS 
analysis. Based on the above considerations, we there- 
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FIG. 5: Temperature dependence of the first cumulant for the 
nearest neighbor single scattering path in Cu, with and with- 
out the perpendicular correction from Eq. [71 and obtained 
either from the minimization of the free energy (FE) or from 
the correlated Einstein model (EM). Both experimental dif- 
ference values witbi^i and withoul*^ the perpendicular motion 
correction were shifted to match the LDA (I) results at K. 

fore propose two alternative prescriptions to stabilize our 
DW factor calculations: 

I. Our first prescription is based on DFT calculations 
using the LDA exchange-correlation functional at the ca/- 
cfiZaierf equilibrium lattice constants a(T) at a given tem- 
perature. Note however that the errors in the LDA esti- 
mates of the lattice constant are often larger than those 
obtained in fits to XAFS experiment. 

IL Our second prescription is based on DFT calcula- 
tions using a modified GGA exchange-correlation func- 
tional termed hGGA (with half-LDA and half-GGA) at 
the experimentally determined lattice constant a(T) at a 
given temperature. As described below, this functional 
is constructed on the assumption that the "true" func- 
tional lies somewhere between pure LDA and GGA. This 
second prescription may be useful, for example, during 
fits of XAFS data to experiment, during which the inter- 
atomic distance is refined. 

Clearly, the use of experimental structural parameters 
limits prescription II, since it requires the knowledge of 
the crystal structure at each of the temperatures of in- 
terest. Such information is only available a priori for a 
handful of systems, although it could be introduced as 
part of the fit procedure. 



B. Exchange-Correlation Functionals 

In the course of this work, we investigated a num- 
ber of exchange-correlation functionals. Generally, the 
exchange-correlation functional is attractive and hence 
strongly affects the overall strength and curvature of the 
interatomic potential. On the other hand it is well known 
that LDA functionals tend to overbind, yielding lattice 



TABLE 1: Born- von Karman parameters DYj (N/m) from 
neutron scattering compared with ab initio calculations from 
this work. 
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constants smaller than experiment typically by about 
1%. In contrast, GGA functionals tend to underbind— 
by about the same amount. These errors are confirmed 
by our calculations, which show that for Cu the LDA 
yields a lattice constant of 3.57 A at K and 3.58 A at 
298 K, while the GGA yields 3.69 A and 3.70 A respec- 
tively, experiment being 3.61 A. Moreover, the effect of 
the functionals on the phonon structure is even larger. 
For example, Narasimhan and de Gironcoli^- show that 
the thermal expansion is about 10% high with LDA and 
10% low with GGA. 

Although significant effort has been put into so-called 
meta-GGA functionals2ii2^ that address these issues, 
they have not yet been widely implemented. Therefore, 
to be consistent with various results^^ and to preserve 
the advantages of the LDA and GGA functionals we have 
devised a modified functional termed hGGA which is a 
mixture of 50% LDA and 50% GGA, i.e., with a 50% re- 
duction in both the additional exchange and correlation 
terms in the GGA. The motivation for the 50-50 mixture 
stems from the observation that the experimental values 
for many quantities are roughly bracketed by the LDA 
and GGA predictions. To simplify the development, we 
chose the closely related Perdew-Wang 92^^ (LDA) and 
Perdew-Burke-ErnzerhoP^ (GGA) functionals. For this 
case the equal parts mixing can be achieved with two 
simple changes: First, the n parameter in the exchange 
energy term in PBE is reduced by half. This change 
preserves all the conditions on which PBE was founded, 
except the Lieb-Oxford bound. Second, the gradient con- 
tribution H to the correlation energy is also reduced by 
half. Similar modified functionals for solids have been 
proposed by Perdew et al,"^^ suggesting that modifica- 
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TABLE II: Debye- Waller factors al{T) (in lO'^ A^) for the 
single scattering path to the n"" shell of some fee lattice met- 
als. The experimental difference values were shifted to match 
the LDA (I) at 80 K and the experimental error (in parenthe- 
ses) indicates the error in the least significant digits. 
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tions similar to the hGGA may be more generally appli- 
cable. Fig. [T] shows the average frequency obtained with 
the hGGA functional and confirms that this yields the 
desired behavior. 



C. Dynamical Matrix 

The key physical quantity needed in calculations of 
the Debye- Waller factors is the dynamical matrix D. 
With modern electronic structure codes this matrix of 
force constants can be calculated with sufficient accu- 
racy from first principles both for periodic and molecular 
systemsi^i^ In this paper we restrict our attention to 
periodic systems which can be treated, for example, us- 
ing the methodology implemented in the ABINIT— code, 
as described in detail in Ref. Briefly, the reciprocal 
space dynamical matrix 

D,o.rp (q) ^Y.D,^o.,,'vpe^''<^^'^'-'^'°) (15) 

is calculated in a 4 x 4 x 4 grid of q- vectors. This grid 
is interpolated inside the Brillouin zone and the real- 
space force constants arc obtained by means of an inverse 
Fourier transform. We find that such an interpolated grid 
gives well converged real-space force constants up to the 
fourth or fifth shell. The neglect of the shells beyond 
that introduces an error much smaller than other sources 
of error in the method. Finally, since the calculation of 
the DW factors uses clusters that typically include about 



TABLE III: Debye- Waller factors cr^(r) (in 10^^ A^) for the 
single scattering path for the n"* shell of some diamond lat- 
tice semiconductors. The experimental difference values were 
shifted to match the LDA (I) at 80 K and the experimental er- 
ror (in parentheses) indicates the error in the least significant 
digits. 



n CD LDA(I) GGA hGGA(II) Expt 



Ge 


1 


5.11 


3.42 


3.98 


3.76 


3.5(1) 


[29] 


295 K 


2 


7.43 


10.38 


11.91 


9.70 


9.6(8) 




3 


7.64 


13.09 


15.03 


11.84 






GaAs 


1 


5.17 


3.97 


4.59 


3.86 


4.2(1) 


[29] 


295 K 


2 


7.75 


12.70 


14.28 


12.01 


11.7(14) 




(Ga Edge) 


3 


7.69 


14.91 


16.29 


14.01 






GaAs 


1 


5.15 


3.96 


4.59 


3.86 


4.2(1) 


[29] 


295 K 


2 


7.20 


10.80 


11.69 


10.19 


9.6(11) 




(As Edge) 


3 


7.68 


14.83 


16.66 


14.00 







20 shells, the full force constant matrix for these clusters 
must be built by replicating the 3x3 Djia.j'i'p blocks 
obtained for each jl,j'l' pair. 

D. Lattice and Force Constants 

The temperature-dependent lattice constant a{T) is 
obtained by minimizing F{a,T) in Eq. ([6]) with respect 
to a at a given temperature T. Within the electronic 
structure code ABINIT, the total VDOS Pa{oj) is calcu- 
lated with histogram sampling in q-space. However, we 
find it more convenient here to use a Lanczos algorithm 
in real space, similar to the approach used for the MSRD. 
This can be done by modifying the initial normalized dis- 
placement state ]0) in Eq. fO]) to that for a single atomic 
displacement, rather than the displacement along a given 
MS path. If more than one atom is present in the unit 
cell the contributions from each atom must be calculated 
and added. Similarly for anisotropic systems one must 
trace over three orthogonal initial displacements. Fig. [2] 
shows a typical VDOS generated using the Lanczos al- 
gorithm. We find the free energies calculated with this 
approach deviate from the q-space histogram method by 
less then 2 meV, i.e., to within 1%. 

To minimize F(a, T) efficiently we proceed as follows: 
First, the lattice constant is optimized with respect to 
the internal energy E{a) and a potential energy surface 
(PES) for the cell expansion is built around the mini- 
mum. Second, the ab initio force constants are computed 
at each point of the PES to obtain the vibrational com- 
ponent of F{a, T). Since this is the most time-consuming 
part of the calculation, we have taken advantage of the 
approximately linear behavior for small variations as il- 
lustrated in Fig. [TJ Then, each element of the force 
constants matrix is interpolated according to 

Djla,j'l'l3 — Ajia^j'1'13 + Bjia^j'l'p Aa (16) 

from just two ah initio force constant calculations with 
slightly different lattice parameters. This interpolation 
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FIG. 6: Temperature dependence of the perpendicular com- 
ponent a\ of the Debye- Waller factor for the nearest neigh- 
bor single scattering path in Cu. For comparison we also plot 
2.5 to show the correlation a\ « 7±(J^, together with the 
values extracted from experiment.— 



scheme allows us to reduce the computational cost of a 
typical calculation by a factor of 2/3, while introducing 
an error of less than 2% in the average frequencies. Once 
the values of F(a, T) on the PES are obtained, we de- 
termine the minimum a{T) by fitting -F(a, T) to a Morse 
potential 



IV. RESULTS 
A. Born-von Karman parameters 

Phonon dispersion curves are often parametrized in 
terms of so-called Born-von Karman (BvK) coupling con- 
stants. These parameters are essentially the Cartesian 
elements of the real space dynamical matrix defined in 
Eq. pil)) . The main difference between the Born-von 
Karman parameters and force constants obtained within 
the quasi-harmonic approximation is that the former are 
tabulated at specific temperatures while the temperature 
dependence of the quasi-harmonic model arises implicitly 
from the dependence of the lattice parameters on ther- 
mal expansion. The dominant BvK coupling constants 
(up to the second neighbor) are presented in Table HI 

We find that both the LDA with prescription I and the 
hGGA with prescription II generally give force constants 
that are within a few percent of experiment. Typically 
the LDA force constants with prescription I are slightly 
higher than those from the hGGA with prescription II. 
Also, note that the transverse components of the BvK pa- 
rameters tend to be overestimated. We have also consid- 
ered the pure PBE GGA functionals, but find that they 
produce force constants that are significantly weaker due 
to their larger equilibrium lattice constants (Fig. [T]). 



B. Mean-square Relative Displacements 



F(a,T) -i?o 



-2p(a-a(T)) _ l3(a-a{T)) 



(17) 



We have estimated that the numerical error in this mini- 
mization is of order 5 x 10^'' A or less by fitting only the 
internal energy component E{a) and comparing with the 
minima obtained using conjugate gradient optimization. 



E. Computational Details 

All the ABINIT calculations reported here use 
TrouUier-Martins scheme — Fritz-Haber-Institut pseu- 
dopotentials. We found that an 8 x 8 x 8 Monkhorst- 
Pack fc-point grid and an energy cutoff of 60 au (12 au for 
Ge) were sufficient to achieve convergence with respect to 
the DW factors. In all cases where the geometries were 
varied, an energy cutoff smearing of 5% was included 
to avoid problems induced by the change in the num- 
ber of plane wave basis sets. For metallic systems, the 
occupation numbers were smeared with the Methfessel 
and Paxton^^ scheme with broadening parameter 0.025. 
Results are presented for LDA (Perdew-Wang 92^^) and 
GGA (Perdew-Burke-Ernzerhof^^) functionals, as well as 
for our mixed hGGA functional. 



Calculations of the MSRD for the dominant first near 
neighbor path for fee Cu are shown in Fig. [31 and de- 
tailed results for various scattering paths are presented 
in Table [Til Both of our prescriptions I and II yield re- 
sults in good agreement with experiment. For Cu even 
the correlated Debye model is quite accurate. Note also 
a slight deviation from linearity in temperature T due 
to the variation in the dynamical matrix with tempera- 
ture is visible both in the experimental curve and in the 
calculation using prescription I. 

Similarly, calculations of the MSRD for the first neigh- 
bor path in Ge are shown in Fig. [4l and detailed re- 
sults for various scattering paths are given in Table [Till 
Again, both of our prescriptions yield results in good 
agreement with experiment, with the LDA prescription 
being slightly better. For this case, however, the corre- 
lated Debye model is significantly in error; this is not 
unexpected given the strong anisotropy of the diamond 
lattice. Tables [iTl and IIIII also include similar results for 
Ag, Pt and GaAs. 



C. Thermal Expansion 

The thermal expansion can now be calculated in two 
ways. First, by minimizing the free energy of the sys- 
tem using Eq. ^ one can obtain the overall thermal 
expansion corresponding to the expansion of the lattice 
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T(K) 



FIG. 7: Temperature dependence of the third cumulant for 
the nearest neighbor single scattering path in Cu. The exper- 
imental difference values^^ were shifted to match the LDA (I) 
results at K. 



low temperature. We also show that for fee structures the 
correction due to perpendicular motion is smaller than 
the crystallographic contribution by a factor of 7_l/67, 
which for Cu is about 20%. To illustrate this correlation, 
Fig. [6] shows the perpendicular motion contribution 
calculated both by the Lanczos procedure and with a 
constant correlation factor 7^ = 2.5. 

We have carried out similar calculations of cr^ for the 
case of diamond lattices. Due to the strongly directional 
bonding in diamond structures, and non-negligible bond 
bending forces, the calculations are more complicated 
than for fee materials. Our ab initio calculations using 
the LDA with prescription I yield a ratio that varies from 
3.4 to 7.2 between and 600 K, in reasonable agreement 
with experiment where 7_l varies between 3.5±0.6 and 
6.5±0.5 in the same range.— In contrast our single near 
neighbor spring model (Appendix) gives a smaller high 
temperature value 7j_ = 3.5 and the addition of a sin- 
gle bond-bending parameter does not improve the agree- 
ment. 



constant a{T). For monoatomic systems the thermal ex- 
pansion of any MS path is simply proportional to the 
lattice constant. More generally, the expansion is MS 
path dependent, and can be estimated using the corre- 
lated Einstein model of III CI and the Appendix. From 
Eq. and the Einstein model Griineisen parameter 
7 = —k^R/k, this model predicts that the first cumu- 
lant a^-^"^ has a temperature dependence proportional to 



3777 
R 



^(1) ^ 



(18) 



As shown in Fig. [5] (dashed and dotted curves), this corre- 
lated Einstein model estimate for the thermal expansion 
agrees well with that obtained from minimizing the free 
energy of the system and with experimental crystallo- 
graphic data. 



D. Perpendicular Motion Contributions 



E. Third Cumulant 

As for the first cumulant, the third cumulant can be 
estimated from the correlated Einstein model, and the 
relation 



= ^ 



2- - ^ 



(19) 



Again an additional scaling factor rj is needed to correct 
the original Einstein model relations when cr'-'^' and cr^ 
are replaced by the full results from our LDA calcula- 
tions. Also the presence of this factor rj gives another 
correction to the relation a^^V^/cr'-^-' « 2 given by clas- 
sical models^ or the correlated Einstein model at high 
temperatures. 



Fig.[n]also shows the first cumulant for Cu obtained by 
adding the crystallographic component and the 

correction due to perpendicular motion Act^^'' from Eq. 
([7]). As observed by Fornasini et aL^^ the mean square 
perpendicular motion (MSPD) is correlated with cr^, i.e., 
CT^ — j±cr'^, with an observed proportionality constant 
for Cu 7j_ « 2.5 ± O.SM The MSPD can be calculated 
using our Lanczos procedure with an appropriately mod- 
ified seed state |0) for perpendicular vibrations. This 
yields a ratio for 7_l that varies from 2.17 to 2.36 be- 
tween and 500 K, respectively, for Cu. Moreover, as 
shown in the Appendix, this ratio can also be estimated 
using a correlated Einstein model for fee structures, and 
we derive a value of 2.5 at high temperatures. The cor- 
related Einstein model also predicts that 7j^ is weakly 
temperature-dependent, reducing to about \/5 ~ 2.24 at 



F. Crystallographic Debye- Waller Factors 

Finally, we present results for the x-ray and neu- 
tron crystallographic Debye Waller factors W{T) — 
(l/2)fc^M^(T), where the mean-square displacements 
M^(T) = ((u • R)^) are given by Eq. with /Oq(w) given 
by the total vibrational density of states per site, as cal- 
culated by our Lanczos algorithm with an appropriate 
seed state. For this case good agreement is obtained 
for both of our DFT prescriptions at low temperature, 
although the errors become more significant at higher 
temperatures. Also, we find that the convergence of the 
Lanczos algorithm is slower than for the path dependent 
Debye- Waller factors, requiring approximately 16 itera- 
tions to achieve convergence to 1%. 
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n 

50 100 150 200 250 300 350 400 450 500 
T(K) 

FIG. 8: Temperature dependence of the crystallographic 
Debye- Waller factor for the nearest neighbor single scatter- 
ing path in Cu, and compared to experimental values^. 



V. DISCUSSION AND CONCLUSIONS 



We have developed a first principles approach for cal- 
culations of the Debye- Waller factors in various x-ray 
spectroscopies, based on DFT calculations of the dynam- 
ical matrix and phonon spectra for a given system. We 
find that the results depend strongly on the choice of 
exchange-correlation potential in the DFT, but we have 
developed two prescriptions that yield stable results, one 
based on the LDA and one based on a modified GGA 
termed hGGA. Calculations for the crystalline systems 
presented here show that our LDA prescription yields 
good agreement with experiment for all quantities, typ- 
ically within about ± 10%. Second, if the lattice con- 
stant is known a priori, our hGGA prescription also pro- 
vides an accurate procedure to estimate the MSRD. An- 
harmonic corrections and estimates of the contribution 
from perpendicular vibrations are estimated using a cor- 
related Einstein model. For these anharmonic quantities, 
however, we have found that the comparative softness 
of the lattice dynamics with the GGA and hGGA func- 
tional leads to results which are somewhat less accu- 
rate than those for the LDA. Finally we have also calcu- 
lated the crystallographic Debye Waller factors. Our ap- 
proach also yields good results for calculations of DW fac- 
tors in anisotropic systems, as illustrated for for Ge and 
GaAs. All of these results demonstrate that the prescrip- 
tions developed herein can yield quantitative estimates of 
Debye- Waller factors including anharmonic effects in var- 
ious crystalline systems, and generally improve on phe- 
nomenological models. Extensions to molecular systems 
are in progress. 
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APPENDIX 

In this Appendix we briefly discuss the correlated Ein- 
stein model used in estimating anharmonic contributions 
to the DW factors. The model is illustrated with an 
application to the correlated Einstein model for calcu- 
lating the mean-square radial displacement (MSRD) 
and mean square perpendicular displacement (MSPD) 
ai = (|Auj.p). ^ 

The construction of Einstein models is not unique in 
that different physical quantities reflect different aver- 
ages over the VDOS. For example, the theoretical MSRD 
given by Eq. ^ reflects an average over a thermal weight 
factor varying as l/oj^ at high temperatures. Thus the 
Einstein model parameters in our prescription are con- 
structed to preserve the correct high temperature behav- 
ior of cr^ji!^ The first step in this construction is the cal- 
culation of from the total potential energy for a net 
displacement a; of a path along a particular seed displace- 
ment state |0). Next this value is renormalized to give 
the correct MSRD at high temperatures. Thus we define 

/co = ?7/i(0|D|0) (A.l) 

where o;^ is given from Eq. ([72]) and the factor ?/ = 
l/(ti;~^)(D^, where is the inverse second moment 

of the projected VDOS. The cubic coupling is then 
defined to be consistent with the variation in k given by 
the Griineisen parameter 

and hence must be is similarly renormalized 

fc3 = |^(0|D(i?)|0). (A.3) 

Then the Einstein frequency a;^; in the quasi-harmonic 
approximation is obtained from the relationi^ 

k = ko + Gk^x = [i<^\. (A. 4) 

where /x is the reduced mass. For Cu using the LDA 
(I) prescription for the dynamical matrix, this procedure 
yields r] = 0.73, fco = 54.7 N/m, fcg = -48.4 N/mA, and 
k = 51.1 N/m. 

The scaling factor rj thus forces the relation a'^(T) 
(J%{T) = kBT/^oj"^ at high temperature, where (j'^{T) in 
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the Einstein model is 



; coth 



(^) 



(A.5) 



and the zero-point value CTq = cr^(O) = hujE/^k. 

Then, from ct^ and relations between the 
cumulants^ii^ one can obtain MS path-dependent 
estimates for cr^^^ and a^^\ When these relations are 
expressed in terms of the calculated (or experimental) 
values of the cumulants, we have found it necessary to 
include multiplicative factors of rj compared to the pure 
Einstein model expressions, ^"^^ to obtain quantitative 
agreement with experimental results e.g., as shown in 
Figs. El and [3 

As a second example, we construct such a model for 
monoatomic fee Cu starting from an anharmonic pair po- 
tential. That is, we will assume that the lattice dynamics 
can be described by an anharmonic pair potential Vb be- 
tween near-neighbor bonds of the form 



Vo{x) ^ ^k°x^ + k^x^ 



(A.6) 



Here x is the net displacement x along the bond direction, 
with positive displacements referring to expansion and 
negative to compression.^;-^? 

First consider the potential energy Vn (x) for vibra- 
tional displacement x along the bond {OR) between 
lattice points (0,0,0) and R = i?(0, 1, 1)/V2. The 
net anharmonic potential (x) is then given by Eq. 
(fTSll with a displacement state x\0) defined by uq = 
(x/2)(0,l,l)/\/2, and = (-x/2)(0, 1, 1)/V2. Then 
constructing the dynamical matrix using Eq. ((X6ll with 
small displacements, we find a net spring constant fco — 
77(0|L»|0) = {5r]/2)k°, in agreement with Ref. 0. This 
result can alternatively be obtained by summing the 23 
pair potentials between the shared bond (OR), the 11 
nearest neighbor bonds to the origin and the 11 other 
nearest neighbor bonds to R, giving V||(a;) — Vo{x) + 
2Vo{-x/2) + 8Vb(a;/4) + 8Voi-x/A) + 4Vb(0). Similarly 
we find that the anharmonic coupling [cf. Eq. (|A.3p ] is 
fca = (377/4) fc!^, so that 



(A.7) 



where we have again included a factor rj so that the Ein- 
stein model for agrees with the expression from the 
inverse second moment of the VDOS. 

In a similar way, we can develop a correlated Ein- 
stein model to describe the perpendicular vibrations. For 
this case we consider a vibrational displacement Au± 
of length y — |Auj^| perpendicular to the bond be- 
tween {0,R). Thus we set uq = (?//2)(0, 1, -1)/V2 and 
uji = (— y/2)(0, 1, — 1)/-\/2, where R is the nearest neigh- 
bor distance. The net potential V±{y) is again obtained 



by summing the 23 pair potentials between the shared 
bond (OR), the 11 nearest neighbor bonds to the ori- 
gin and 11 others to _R, similar to the calculation above 
for vibrations along the bond. For this case, two bonds 
are stretched by y/2, two contracted by y/2, three un- 
changed, eight stretched by y/A, and eight contracted by 
2//4, yielding a net sum V±{y) = 2Vo{y/2) + 2Vo{^y/2) + 
8Vb(y/4) + 8Vb(-?//4) + 3Vb(0), and hence 



(A., 



Note that by symmetry, the net cubic anharmonic con- 
tribution vanishes. Thus the effective spring constant for 
the MSPD is fc = 2rik'^ and predicted to be insensitive to 
thermal expansion. The correlated Einstein model Vj_(z) 
is clearly the same for the MSPD along the 2;-axis, 

With these results we can show that the MSPD for 
the first neighbor path in fee materials a'j_ is correlated 



with (T^ = (\Ai 



Both the MSPD and MSRD in the 



Einstein model are given by Eq. lA.Si with their respec- 
tive Einstein frequencies uje — {k/ny/^. For the total 
contribution from perpendicular vibrations, one has to 
multiply by two to get the net cr^ from both independent 
axes. At high temperatures for example, we obtain the 
MSPD al = 2kBT/2rik°. This is higher than the MSRD 
cr2 = {2/5)kBT/r]k° by a factor of 7_l = 5/2. The model 
also predicts a weakly temperature dependent ratio 

^ (T\ = din = o ^ECothi/3ncjj,/2) 

' a2(T) iu^coth{fihLUE/2)- ^ ' 

This ratio varies between -y/S — 2.236 and 2.5 with in- 
creasing temperature. Thus the ratio j± obtained with 
the correlated Einstein model, for the fee lattice depends 
only on geometry and describes the anisotropy of the vi- 
brational ellipsoid in monoatomic fee structures reason- 
ably well. 

Because of the above relation between ct^ and cr^ in 
the Einstein model, the perpendicular motion correction 
can be related to the contribution to lattice expansion 
from anharmonicity. Thus from Eq. ([7]) and Eq. ^TE\\ . we 
find that 



7_l/6777. 



(A.IO) 



For fee Cu this ratio predicts a correction to the first cu- 
mulant a^^^ from perpendicular motion of about 25%. In- 
deed, this shift is comparable to the observed differences 
in the thermal expansion with and without the perpen- 
dicular motion correction observed in Fig.[5l Thus for the 
dominant near neighbor bonds, the correlated Einstein 
model predicts a comparatively small but non-negligible 
effect of perpendicular motion on EXAFS distance deter- 
minations. 
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